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Approximation scheme for master equations: variational approach to multivariate case 
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We study an approximation scheme based on a second quantization method for a chemical master 
equation. Small systems, such as cells, could not be studied by the traditional rate equation approach 
because fluctuation effects are very large in such small systems. Although a Fokker-Planck equation 
obtained by the system size expansion includes the fluctuation effects, it needs large computational 
costs for complicated chemical reaction systems. In addition, discrete characteristics of the original 
master equation are neglected in the system size expansion scheme. It has been shown that the usage 
of the second quantization description and a variational method achieves tremendous reduction in the 
dimensionality of the master equation approximately, without loss of the discrete characteristics. We 
here propose a new scheme for the choice of variational functions, which is applicable to multivariate 
cases. It is revealed that the new scheme gives better numerical results than old ones and the 
computational cost increases only slightly. 

PACS numbers: 82.20.-w, 02.50.Ey, 05.40.-a, 87.10.-e 



I. INTRODUCTION 

Time evolution of chemical reaction systems can be 
studied by a traditional rate equation approach. In the 
rate equation approach, all fluctuation effects are ne- 
glected. In order to treat the chemical reaction systems 
more appropriately, it is possible to write discrete mas- 
ter equations which include all information of the fluc- 
tuation effects. In recent years, it has been necessary to 
study a stochastic process in a system with small size. 
Here, the word 'small size' means that the number of 
molecules contained in the system is very small. In such 
small size system, the fluctuation effects are large and 
cannot be ignored. For example, some of biochemical 
systems, such as gene regulatory networks, contain small 
number of regulatory molecules [ij, [2J . It is expected that 
the fluctuation effects and the discrete effects of the orig- 
inal master equation play an important role in such small 
size systems. 

In general, we cannot solve the chemical master equa- 
tions analytically. While one might think that direct nu- 
merical evaluation of the chemical master equation is the 
best way in order to investigate such small systems, it 
would be also impossible. The master equation can be 
understood as a large set of coupled ordinary differential 
equations, and the number of them often becomes very 
large; it increases exponentially with the number of chem- 



ical substances. For example, we have 100 
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pled differential equations when we consider a small size 
system of only five chemical substances with the num- 
ber of each chemical substances varying between and 
99. It is possible to approximate the chemical master 
equations by a Fokker-Planck equation which is obtained 
by the system size expansion method 0, HJ- However, 
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the Fokker-Planck equation has the similar problem of 
the curse of the dimensionality, i.e., the computational 
costs increases exponentially with the number of chem- 
ical substances. In addition, the system size expansion 
method neglects the discrete characteristics of the orig- 
inal chemical master equation, and the system size ex- 
pansion method becomes a good approximation only if 
we consider a chemical system with large size. 

A chemical reaction system consists of Markov jump 
processes in continuous time. Such dynamics can be sim- 
ulated (exactly) using several Monte Cairo methods such 
as the Gillespie algorithm [§]. In the Gillespie algorithm, 
the lapse time to the next event is determined by expo- 
nentially distributed random numbers, and we determine 
which event occurs in proportion to the rate of the event. 
While the numerical simulations are available for study- 
ing complicated stochastic systems, these methods could 
sometimes become time-consuming when there are many 
chemical reactions, because averaging procedures of a lot 
of Monte Carlo samples are needed in order to obtain 
statistics values. Hence, it is important to study analyt- 
ical methods and approximation schemes applicable to 
small systems. 

The field theoretic approach (or second quantization 
description) is one of the candidates for such analytical 
treatments. The analogy of the master equation to quan- 
tum descriptions has been introduced by Doi, and several 
authors have developed the formalism p, 7, 8]. The field 
theoretic approach has revealed the anomalous kinetics in 
reaction-diffusion systems incorporating the renormaliza- 
tion group method 0, [l(| , and the analytical scheme has 
been applied to various phenomena [ll|, [H, HH, EH • 
In addition, the variational method based on the second 
quantization description developed by Sasai and Wolynes 
[2| is hopeful in order to investigate the fluctuation and 
discrete effects in the small size systems. The variational 
method is an approximation, and then we can avoid the 
curse of dimensionality. In addition, the method does not 
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ignore the discrete characteristics of the original chemical 
master equations. 

The variational method has been already applied to 
stochastic gene regulatory networks, and it has given 
qualitatively good results [3, HB, 113] • However, the origi- 
nal formulation of the variational method has been based 
on the Poisson ansatz, in which we assume that the prob- 
ability distribution of the system is described by the Pois- 
son distribution. Hence, it is necessary to use a Hartree 
approximation, and then only restricted fluctuation ef- 
fects can be included |17| . In order to improve the 
variational method, "a superposition ansatz" has been 
proposed in ref. [l8|. The new ansatz is based on the 
coherent state, and gives a quantitatively good result for 
master equations describing the gene regulatory network. 

The main purpose of the present paper is to investi- 
gate the advantage of the new variational ansatz and its 
application limitation. In ref. the validity of the 

superposition ansatz has been confirmed for a one-body 
problem in which the Hartree approximation is adequate. 
However, the applicability of the superposition ansatz 
for multivariate cases has not been studied. For practi- 
cal purposes, the multivariate case is important for real 
chemical reactions in cell. Hence, in the present paper, 
we mainly focus on the applicability of the new ansatz to 
multivariate problems. We first consider a simple uni- 
variate problem in order to investigate zero-boundary 
effects: the reaction system stops when the number of 
molecules becomes zero. For the univariate problem, we 
clarify that the survival probability is easily calculated 
from the Laplace transformation of the variational func- 
tion. Next, we consider a bivariate problem which cannot 
be treated by the original Poisson ansatz adequately. For 
the bivariate problem, we will show that the Hartree ap- 
proximation does not give quantitatively good results in a 
chemical reaction system with small size. The superposi- 
tion ansatz gives better results than those of the Hartree 
approximation. In addition, it enables us to calculate the 
covariance, which is incomputable by the Hartree approx- 
imation. 

The construction of the present paper is as follows. 
In Sec. II, we will explain the formalism of the second 
quantization method and the variational method. Sec- 
tion III gives numerical results for a univariate problem, 
and we will show an easy calculation scheme for the sur- 
vival probability. In Sec. IV, we treat a bivariate problem 
and study the availability of the superposition ansatz. 
Section V gives concluding remarks. 



II. FORMALISM 

At first, we summarize the second quantization de- 
scription for a master equation and the variational 
method developed in refs 0] and [lil ]. 



A. Second quantization description 

We consider a system with d chemical substances, 
X\ , . . . , Xd, and denote the number of each chemical sub- 
stance Xi as rii. Each reaction has the following generic 
structure: 

3 3 

and hence the probability of a transition per unit time is 
proportional to the rate c a and the number of molecules, 
giving a reaction rate of 

where {n-i} = {ni,n2, ■ ■ ■ In this case, the master 

equation for the probability of a configuration {rii} has 
the form: 

j t p ({ni}) = £^ a (K}) p (K}) - i? a (K})p(K}), 

a 

(3) 

where nf = n, + vf — rjf is the number of molecules prior 
to reaction a that leads to a current number rii. 

In general, it is useful to rewrite the master equa- 
tion ([3]) by creation and annihilation operators and to 
treat it by techniques developed by quantum mechanics. 
First of all, we define a ket vector as the state in 

which the configuration of the system is {rii}. We here 
introduce the bosonic operator algebra 

[ai,aj] = Sij, [a,i,aj] = [at, at] = 0, (4) 

where a\ is the creation operator for chemical substance 
Xi, and ai the annihilation operator for Xi. The creation 
and annihilation operators act on the ket vector = 
\ni, . . . , rid) as follows: 

at|m, ...,ni,...,n d ) = K, . . . , + 1, . . -,n d ), (5) 
rai, ...,m,...,n d ) = rn\n\, . ..,rii - 1,.. . ,n d ), (6) 
ata;|ni, . . . , n*, . . . ,n d ) = rii\ni, . . . ,n i5 . . . ,n d ). (7) 

Therefore a state with the configuration {ni, . . . ,n d } is 
obtained from the empty vacuum state |0) as 

d 

IW>-II( a I) n! l )' ( g ) 

where the empty vacuum state |0) is defined as Oj|0) = 
for all i. 

Using the time-dependent probability distribution of 
the original discrete master equation, P({rii}), we intro- 
duce the ket vector of the system as 

|*> = E p (K»IW>- ( 9 ) 
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The original master equation is then translated into the 
following equation 



0_ 

Of 



(10) 



where f2 is the time evolution operator. The time evo- 
lution operator O is constructed so that it becomes con- 
sistent with the original master equation. We will show 
concrete examples of the time evolution operator in 
Sees. Ill and IV. 

We remark on the computational method for the calcu- 
lation of averages with respect to the original probability 
distribution P({m}). The average is taken by the usage 
of the projection state (V\ = (0| IliLi exp(ai). The pro- 
jection state satisfies (V^) = 1 and a statistical average 
of an observable A is obtained by 



(A) = ]T A{{n t })P{{ ni }) = (^({otoi})!*) 
{«>} 



(11) 



Using the above theoretical scheme, we can calculate 
various quantities for the original chemical master equa- 
tion rigorously. However, the evaluation of the time de- 
pendent ket vector \ty) is a computationally difficult task, 
and in general it is impossible to obtain the exact solution 
of the ket vector. We, therefore, use an approximation 
method in order to obtain the time-evolution of the ket 
vector. 



B. Variational method 

In order to reduce the dimensionality of the problem, a 
variational method developed by Eyink 0, [2(| is avail- 
able. We here briefly review the variational method com- 
bined with the second quantization description Q ■ 

When we define an effective action T as 



r = y dt<$|(fit-n)|¥>, 



(12) 



Eq. (fT0|) is equivalent to the functional variation = 
0. Because of the non-Hcrmitian property, it is not al- 
ways true that the left eigenvectors and right eigenvectors 
are the same. Hence, we assume two variational functions 
for the bra and ket states, respectively: The ket state \ty) 
(or the bra state (<E>|) is parametrized by a R (or a L ), and 
where a R and a L are vectors with K components; 



(13) 
(14) 



A set of finite dimensional equations for parameters a R 
and a L is obtained by the functional variation procedure. 
Note that we set $(a L = 0) to be consistent with the 
probabilistic interpretation, so that 



($(a L =0)|*(a R )) = 1. 



(15) 



We, therefore, obtain the following equation which stems 
from an extremum of the action 



K 

E 



dak, 



\ da? 
~df 



da? 



<9$ 
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•>] 







= 



for m = 1, 2, • • • , K. 

(16) 

Using this variational scheme, we obtain a set of time- 
evolution equations for the time-dependent parameters 
a R , and the equations can be solved numerically. The 
only remaining task is to give an explicit ansatz for ($| 
and I*). 



C. Poisson ansatz 

The Poisson ansatz has been introduced in ref. 0. In 
the Poisson ansatz, we assume that the probability of the 
number of each molecule, rii, obeys the Poisson distribu- 
tion. To be more precise, the probability distribution of 
the original discrete master equation assumes to have a 
product form 



P({m}h 



rid\ 



(17) 



where fii is the mean value of the number of molecule Xi. 
This means that we need a Hartree approximation 

i*) = \ipi)®---®m, (is) 

and each one-body ket vector \ipi) is assumed as 



exp 



Mi (a! 



1 



10,), 



(19) 



where |0j) is the vacuum state for a chemical substance 
Xi. In this case, we have d time-dependent parameters 
q r = {/ii, . . . , fid}, and by using an adequate bra ansatz, 
d coupled differential equations for the time-dependent 
parameters {ni} are obtained. 

We note that the Poisson ansatz includes only re- 
stricted fluctuation effects because the mean value and 
the variance of the Poisson distribution are the same. In 
addition, we cannot calculate any correlations between 
chemical substances because the Poisson ansatz needs the 
Hartree approximation. 



D. Superposition ansatz 

We here use the analogy between the Poisson ansatz 
and the coherent state of quantum mechanics. A new 
ansatz, "superposition ansatz," is constructed by the su- 
perposition of the coherent states as follows: 



poo d 

|*) = / dxh({xi};{^})T[exp[xj(a\ 
Jo 



1)]|0), (20) 
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where J °° dx = / °° dx\ . . . J °° dxd, and h({xi}; {/x}) is 
the variational function. Here, {fi} is a set of the varia- 
tional parameters which specifies the variational function 
h({xi}; {fi}). Note that it is possible to use a continu- 
ous probability distribution as the variational function 
h without loss of the discrete characteristics of the orig- 
inal problem, via the usage of the coherent states. In 
general, the continuous property sometimes makes ana- 
lytical treatments easier than discrete cases. In addition, 
the superposition ansatz enables us to use a multivariate 
variational function h({xi}; {fJ-}) so that the correlations 
between the chemical substances can be adequately cal- 
culated by the new ansatz. 



III. EXAMPLE I: UNIVARIATE CASE 

First, we study a univariate case. In this simple ex- 
ample, we mainly focus on the zero-boundary effect: the 
number of molecules cannot be below zero. When we 
use the Fokker-Planck equation, the additional boundary 
condition is needed in order to include the zero-boundary 
effect. In contrast, the variational method does not need 
to consider the zero-boundary, because the superposition 
ansatz is based on the Poisson distributions (coherent 
states) and hence the discrete characteristics is not ne- 
glected. In addition, we will show that the survival prob- 
ability is easily calculated by the Laplace transformation 
of the variational function in the superposition ansatz. 



Model 



We consider the following reaction system: 



X 2X, 
X^0. 



The master equation is given by 



(21) 



^^n, =c\ [{n x - l)P naj -i - n x P n J 

+ c 2 [{n x + l)P n=c+ i - n x P n J 



(22) 



where n x is the number of molecules X. Using the second 
quantization description, we finally obtain the following 
time evolution operator: 



fi = ci(alala x - a x a x ) + c 2 (a x - a x a x ). 



(23) 



It is easy to check the above time evolution operator ad- 
equately recovers the original master equation (f2"2")l . 

We here use the gamma distribution as the continu- 
ous variational function h({x.i}, {fJ,}) in Eq. (|2"01) . The 
gamma distribution is defined as 



F(x;k,9) = x 



fe _! exp(-z/fl) 
T{k)9 k ' 



(24) 



and the mean and the variance are given by kd and k9 2 , 
respectively. This choice of the variational function in- 
dicates that the variational parameters {/x} are given by 
{^1,^2} = {k,0}. The first and second moments of the 
gamma distribution are calculated as 



(x) = k6, 

(x 2 ) = k6 2 + k 2 9 2 , 



(25) 
(26) 



respectively. Note that the above first and second mo- 
ments does not correspond to those of n x . The gamma 
distribution is introduced for the variational function, 
and hence the real probability distribution is obtained by 
the superposition of the Poisson distributions weighted 
by the gamma distribution. The explicit forms of the 
mean and variance of n x will be given in the next sub- 
section. 

By using the following ket ansatz and bra ansatz 



1*)=/ dxF(x-k,9)e^[x(al-l)]\0 x ), (27) 
Jo 

($| = (CMexpj^ + AW^ + A^a,) 2 }, (28) 
the variational parameters are 



a 



= {k,9}, 

= {A«,A( 2 )}. 



(29) 



From Eq. (|16p . we finally obtain the following coupled 
differential equations: 



dk d d6 d 

~dtdk^ X ' + dtd9 
dkd_, 2) Md_ 
dt dk ^ ' + dt 80 



(x) = c x {x) - c 2 (x), 



(30) 



(a; 2 ) = a (2(x 2 ) + 2{x)) - 2c 2 (x 2 ). 

(31) 



Note that the first and second moments, (x) and (a; 2 ) 
depends on the variational parameters k and 9, as de- 
scribed in Eqs. (f2"S"| and (121?!) . and then the derivatives, 
e.g., d(x 2 )/dk, are calculated explicitly. 



B. Numerical result 

We compare numerical results obtained by the super- 
position ansatz with those of the Gillespie algorithm. For 
simplicity, we here set the reaction rates as C\ = 0% = 1. 
The initial parameters are set as fci n i = 30.0 and = 
0.1. When we use the variational method, the only nec- 
essary thing is to evaluate the time evolution of the vari- 
ational parameters k and 9 by using Eqs. (|30|) and (|3ip . 
However, in the Gillespie algorithm, the averaging pro- 
cedure for the Monte Carlo results is needed and we took 
the averages over 10 5 runs. 

Because we set c\ = c 2 , the mean (n x ) does not change 
with time. On the other hands, the variance increases in 
proportional to time. From the variational method, we 
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FIG. 1: (a) The time evolution of the variance (Eq. 
(b) The survival probability (Eq. p6|l '). 



obtain the time evolution of the variational parameters 
k and 9 numerically (from Eqs. ([30]) . (|3"Tj) ). Using these 
variational parameters, the variance of n x is calculated 
as 

Var[n x ] = (n 2 x ) - (n x ) 2 = (x) + (x 2 ) - (x) 2 = k9 + k6 2 , 

(32) 

because 



= VnM dxF(x;k,9) 
Jo 

dxx(l + x)F(x; fc, ( 



{X + (X 



(33) 



In the similar way, it is easy to show the following rela- 
tionship: 

(n x ) = (x). (34) 

Figure Q] shows the numerical results obtained by the 
variational method and the Gillespie algorithm. The 
variance of n x is shown in Fig. [IJa). One can see that 
the results obtained by the variational method and the 
Gillespie algorithm are in good agreement. Note that the 
time evolution of the variance cannot be calculated by 
the Poisson ansatz, because the variance of the Poisson 
ansatz is restricted to the same value as the mean. Next, 



we calculate the survival probability. The survival proba- 
bility is defined as the probability with which the number 
of molecules is not zero. In the variational method, it is 
easy to calculate the survival probability: In a Poisson 
distribution with mean x, the probability with which the 
number of molecules is zero is given by 



. x 
~kl 



(35) 



fc=0 



In the superposition ansatz, the whole probability distri- 
bution is given by the superposition of the Poisson distri- 
bution with mean x weighted by the variational function. 
Hence, by using the Laplace transformation of the varia- 
tional function F(x; k, 9), we obtain the survival proba- 
bility as follows: 



(36) 



P surv = ^ dxF{x; k, 9)e- x = (l + 1 



In Fig. QJb) , we show the numerical results of the survival 
probability. We get a qualitatively good result from the 
variational method; the survival probability shows slow 
decay. In the variational scheme, we construct the proba- 
bility distribution based on the superposition of the Pois- 
son distribution weighted by the gamma distribution, so 
that the probability distribution is not an exact one. By 
using the variational method, a tremendous reduction of 
the dimensionality has been achieved, but there is no free 
lunch here; the price paid is that we cannot obtain an ex- 
act probability distribution. However, unlike the system 
size expansion method, the variational method does not 
neglect the discrete characteristics of the chemical reac- 
tion systems, and hence there is no need to take care of 
the zero-boundary effect as a special boundary condition. 



IV. EXAMPLE II: BIVARIATE CASE 

In the next example, we treat a little complicated reac- 
tion system with two chemical substances. In this case, 
the Hartree approximation is not valid, and we will show 
that the method beyond the Hartree approximation is 
necessary for small reaction systems with large fluctua- 
tion. 



Model 



We consider a reaction system constructed as 






^X, 


X 




Y 










2Y, 



(37) 



6 



The corresponding master equation is 
d 

+ c 2 [{n x + l)(n y - l)P na 

+1,7% —1 W'x'ft'y* n® ,n v \ 
+ c 3 [(% + i)Pn„,n v +l - nyPn*,n y ] 

+ C4 [-Pns.ny-l ~ Pn x ,n v ] ■ (38) 

It is possible to make the corresponding determinis- 
tic rate equations, and the averages of n x and n y ob- 
tained by the deterministic rate equations are (n x ) ra tc = 
c 1 c 3 /(c 1 c 2 +c 2 c 4 ) and (n y )^ atc = (ci+c 4 )/c 3 , respectively. 

The time evolution operator in the second quantized 
description is 

fi =ci(«4 - 1) + (^(o-boJoToj, - a x a x a y a y ) 

+ c 3 (a y - ala y ) + c 4 (a], - 1). (39) 

In the following subsections, wc explain the Hartree 
approximation with two lognormal distributions and a 
multivariate superposition ansatz with a bivariate log- 
normal distribution. 



1. Hartree approximation 

In the Hartree approximation, we set the variational 
function by using two lognormal distributions as follows: 



1*)=/ dx I dyf(x;p x ,a x )f(y;p y ,a y ) 
Jo Jo 

x exp[z(at - l)]exp[y( a t - 1)]|0 X ,0„>, (40) 
($| =(0 X , y \ exp [a x + \ x 1] a x + A<?> {a x f 

+A«a y + A( 2 )K) 2 }, (41) 

where f(x;fi x ,a x ) and f(y]p, y ,a y ) are lognormal distri- 
butions: 



f(x;fx x ,a x ) 



1 



xa x v2n 



exp 



(In a: - p x ) 2 
2ul 



(42) 



and f(y; fi y , a y ) is defined as the similar manner. The 
first and second moments of the lognormal distribution 
are given by 



exp [p x + o- 2 /2] , 
= exp [2p x + 2cr 2 ] . 



„2\ 



(43) 
(44) 



In summary, the variational parameters are as follows: 
a R = {p x ,a x ,p y ,a y }, 



— { A x ^ , Ax \ Xy , Xy }, 



(45) 



and then we obtain the following four coupled differential 
equations from Eq. (I16|) : 



dp x d da x d 

(x) H — - — (x) = ci - c 2 {xy), 



dt dp x 



dt do , 



(46) 



^ X 9 (x 2 ) + ^^(x 2 ) = 2 Cl (x) - 2 C2 (x 2 y), (48) 



cit dp x 
dfiy d 
dt dp y 



(y 2 



dt da x 
da y d 
dt 9cr, 



(y 2 



= 2c 2 (xy 2 ) + 2c 2 (xy) - 2c 3 (y 2 } + c 4 (y). 



(49) 



Since we here use the Hartree approximation, x and y 
do not correlate. Hence, we set (xy) = (x)(y), (x 2 y) = 
(x 2 )(y), and so on. 



2. Bivariate lognormal distribution 

In the reaction system with two chemical substances, 
it is expected that the correlation between n x and n y 
plays an important role, so that we must go beyond the 
Hartree approximation. While the Poisson ansatz needs 
the Hartree approximation, the superposition ansatz en- 
ables us to treat the correlation between n x and n y . In 
order to treat such correlations, we here use the bivariate 
lognormal distribution: 



f(x,y;p x ,p y ,a x ,<7 y ,p) 



2-Kxya x a y \/\ - p 2 



x exp 



1 f 1 



2 Vl-P 2 



In x — p a 



2p 



In x - p x \ / In y - p y 



In V ~ My 



(50) 



The ket ansatz and bra ansatz are set as 

poo poo 

I*) = / dx \ dyf(x,y;fi x ,a x ,fj,y,a y ,p) 
Jo Jo 

x exp[a;(a x - 1)] exp[y(a* - l)]\0 x ,0 y ), (51) 



and 



($1 



={0 XJ W | exp \a x + X x 1] a x + X {2) (a x ) 2 

+X y 1) a y + X y 2) (a y ) 2 + X xy a x a y } 



(52) 
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respectively. By using the bivariate lognormal distribu- 
tion, we can easily obtain the following quantities which 
are necessary to evaluate the time evolution equations: 



(a) 9 



(xy) = exp [pa x a y + p x + a 2 x /2 + p y + o%/2] , 
(x 2 y) = exp [2pa x a y + 2(i x + 2a 2 + p y + a 2 /2] 
(xy 2 ) = exp [2pa x a y + fj, x + a 2 x /2 + 2p y + 2a 2 ] 

The variational parameters are 

a R = {P-x,CTx,P.y,Cry,p}, 

n h _ r A (l) x (2) x (2) , i 

OL — \A X j A x , A y , Ay ,A xy f, 



(53) 
(54) 
(55) 



(56) 



and then we obtain totally five coupled differential equa- 
tions; four equations are equivalent to Eqs. ([4*6")) - ([41?)) . and 
we obtain an additional equation: 



dp x d 



dt dpL x 



(xy) 



da T d 



(xy) 



dp. 



" 9 (xy) 



dt dp, v 



^ — (xy) 
dt day m 



■{xy) 



dt da x 
dp d 
dt dp 1 

= ci(y) + c 2 (x 2 y) - c 2 (xy) - c 2 (xy 2 ) - c 3 (xy) + c 4 (x). 

(57) 

B. Numerical result 

We performed numerical experiments in order to clar- 
ify the advantage and application limitation of the varia- 
tional method. Here, we set c\ = C4 = 1, C2 = 0.01, C3 
i.l and initial parameters are set to p x = p y = log 3, 
a x — a y — 0.1, and for the bivariate lognormal distribu- 
tion, ,o = 0. As in the case of Sec. III. A, the numerical 
results obtained by the Gillespie algorithm are averaged 
over 10 5 runs. 

Figure [2] shows the time evolution of the mean of n x . 
The mean of n x is obtained by 

00 00 />OC poo 

^2 ^2 n x dx dyf(x,y,fx x ,iJ, y ,a x ,ay,p) 

, _n n JO JO 



rij.— rinj—0 



e x x n * e v y Uy 



n x =0 

(x), 



^2 n x dxf(x,y;p x ,p y ,a x ,a y ,p)- 



(58) 



for the multivariate superposition ansatz. For the 
Hartree approximation scheme, we have the same rela- 
tionship; (n x ) = (x). From Fig. [H a )> on e can see that 
the result obtained by the Hartree approximation slightly 
deviates from the results of the Gillespie algorithm. From 
the deterministic rate equations, we predict the mean as 
(n-x)ratc = 5. Hence, the results by the Gillespie algo- 
rithm and the bivariate lognormal distribution suggest 
that the correlation between n x and n y varies the mean 
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FIG. 2: (a) Time evolution of the average of n x . The results 
obtained by the Hartree approximation slightly deviates from 
the results by the Gillespie algorithm, (b) Time evolution 
of the variance of n x . (c) Time evolution of the covariance. 
Here, we set ci = C4 = 1, C2 = 0.01, C3 = 0.1. 



value of n x from that of the deterministic rate equation 
approach. Figure[2jb) shows the variance of n x . Var[n x ] 
is calculated from the same procedure as Eq. ([5!2)) . As 
in the case of the mean, the results obtained by the bi- 
variate lognormal distribution are quantitatively in good 
agreement with the Monte Carlo results. In addition, 
the covariance between n x and n y cannot be calculated 
by the Hartree approximation and, of course, the Poisson 
ansatz. The superposition ansatz enables us to calculate 
the covariance by Cov[n x , n y ] = {xy) — (x) (y). The result 
by the bivariate lognormal distribution is quantitatively 
in good agreement with that of the Gillespie algorithm 
in the present case, as shown in Fig. [!2Jc) . 
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deviation is small compared with the case of the Hartree 
approximation. This means that in such small n x region, 
the higher correlation should be taken into account. We 
expect that the usage of more complicated variational 
function enables us to improve the deviation. 



CONCLUDING REMARKS 
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FIG. 3: Time evolution of the average n x . We here set 

Cl = C 2 = C 3 = 1. 



We showed that the variational method beyond the 
Hartree approximation is useful for calculating physical 
quantities, especially for small systems in which corre- 
lations and fluctuations play important roles. By using 
the variational method, there is no need to take care of 
the zero-boundary, as explained in Sec. III. Of course, it 
is needed to restrict the variational function into some 
specific form, and then it could be difficult to treat the 
correlations and fluctuations correctly in some cases. Fig- 
ure [3] shows the time evolution of the average n x in the 
case with c\ = c% = C3 = 1. The average n x derived 
by the Hartree approximation is largely different from 
the results by the Gillespie algorithm. In addition, the 
results by the bivariate lognormal distribution also devi- 
ates form that of the Gillespie algorithm, although the 



In summary, we clarified that it is possible to apply 
the superposition ansatz to multivariate cases without 
the usage of the Hartree approximation. The superposi- 
tion ansatz enables us to treat various correlations, and 
it gives better numerical results than the Poisson ansatz. 
Such correlations can be included by adding only a few 
differential equations (in our case, only one equation is 
added), and then the computational cost increases very 
little. In addition, we showed that the survival probabil- 
ity is easily calculated via the Laplace transformation. 

In the present paper, we considered a reaction system 
with two chemical species, but in principle, the super- 
position ansatz is applicable to any multivariate case. 
In addition, we believe that the second-quantization 
scheme has a possibility to create more powerful analyti- 
cal method, with the aid of various techniques developed 
in the quantum physics. 

Although we have checked that the multivariate log- 
normal distribution is easy to treat in the variational 
method, it will sometimes be necessary to use more com- 
plicated probability distributions as the variational func- 
tion. As future works, such improvement of the varia- 
tional function is needed. 



[1] M. B. Elowitz, A. J. Levine, E. D. Siggia, and P. S. Swain, 

Science 297, 1183 (2002). 
[2] M. Sasai and P.G. Wolynes, Proc. Natl. Sci. USA 100, 

2374 (2003). 

[3] H. Risken, The Fokker-Planck Equation 2nd edition 
(Springer, Berlin, 1989). 

[4] C.W. Gardiner, Handbook of Stochastic Methods 3rd edi- 
tion (Springer, Berlin, 2004). 

[5] D. T. Gillespie, J. Phys. Chem. 81, 2340 (1977) 

[6] M. Doi, J. Phys. A: Math. Gen. 9, 1465 (1976). 

[7] M. Doi, J. Phys. A: Math. Gen. 9, 1479 (1976). 

[8] L. Peliti, J. Physique 46 1469 (1985). 

[9] D.C. Mattis and M.L. Glasser, Rev. Mod. Phys. 70, 979 
(1998). 

[10] U.C. Tauber, M. Howard, and B.P. Vollmayr-Lee, J. 

Phys. A: Math.Gen. 38, R79 (2005). 
[11] C. Pigorsch and S. Trimper, Phys. Lett. A 300, 221 



(2002). 

[12] R. Dickman and R. Vidigal, Braz. J. Phys. 33, 73 (2003). 
[13] R. Dickman and R. Vidigal, J. Phys. A: Math. Gen. 35, 
7269 (2002). 

[14] J.F. Stilck, R. Dickman, and R. Vidigal, J. Phys. A: 

Math. Gen. 37, 1145 (2004). 
[15] MA. Buice and J.D. Cowan, Phys. Rev. E 75, 051919 

(2007). 

[16] K.Y. Kim and J. Wang, Comp. Biol. 3, 565 (2007). 
[17] K.Y. Kim, D. Lepzelter, and J. Wang, J. Chem. Phys. 

126, 034702 (2007). 
[18] J. Ohkubo, J. Stat. Mech. P09017 (2007). 
[19] G.L. Eyink, Phys. Rev. E, 54 3419 (1996). 
[20] F.J. Alexander and G.L. Eyink, Phys. Rev. Lett. 78, 1 

(1997). 



